library(Zelig)
library(Amelia)
library(texreg)
library(sandwich)
library(dplyr)
options(stringsAsFactors=FALSE)
db = src_sqlite('replication_20160118.sqlite')
bh_1y = tbl(db, sql('SELECT * FROM bh_1y'))
bh_4y = tbl(db, sql('SELECT * FROM bh_4y'))
is = tbl(db, sql('SELECT * FROM investment_shares'))
bh_1y = data.frame(collect(bh_1y))
bh_4y = data.frame(collect(bh_4y))
is = data.frame(collect(is))

# INCOME LEVELS
income = is %>% group_by(income) %>%
         select(matches('orbis|fdi')) %>%
         summarise_each(funs(sum)) %>%
         select(income, orbis_inflow, fdi_in, orbis_outflow, fdi_out)
income = income[c(1, 4:2),]
write.csv(income, 'tables/table_3.csv', row.names=FALSE)

# ESTIMATION FUNCTION
estimate = function(dv, controls, number, data){
    dict = c('bureau'='Bureaucracy', 'corrupt'='Corruption',
             'demoacct'='Accountability', 'ethnic'='Ethnic Conf.',
             'extconf'='External Conf.', 'intconf'='Internal Conf.',
             'invprof'='Investment Profile', 'law'='Law and Order',
             'military'='Military', 'religious'='Religious')
    names(dict) = paste('prs_', names(dict), sep='')
    vars = sort(dict)
    # estimate
    f_apolitical = update(controls, paste(dv, '~ .'))
    f_political = update(f_apolitical, '. ~ . + political')
    models = list(lm(f_apolitical, data))
    for(v in names(vars)){
        data$political = data[, v]
        models = c(models, list(lm(f_political, data)))
    }
    # write to file
    names(models) = c('Baseline', vars)
    fn1 = paste('tables/table_', number, '.txt', sep='')
    fn2 = paste('tables/table_', number + 1, '.txt', sep='')
    idx1 = c(1, 2:6)
    idx2 = c(1, 7:11)
    se = lapply(models, function(x) sqrt(diag(sandwich(x))))
    screenreg(models[idx1], digits=3, omit.coef='iso3c|Intercept|region', stars=NULL,
           override.se=se[idx1], file=fn1)
    screenreg(models[idx2], digits=3, omit.coef='iso3c|Intercept|region', stars=NULL,
           override.se=se[idx2], file=fn2)
}

# BASELINE
estimate(dv = 'unctad_fdi_flow', 
         controls = ~ gdp + growth + trade + inflation + iso3c + year + year2,
         number = 1,
         data = bh_4y)

# STOCK
estimate(dv = 'unctad_fdi_stock', 
         controls = ~ gdp + growth + trade + inflation + iso3c,
         data = bh_4y,
         number = 9)

# PPE
estimate(dv = 'bea_ppe', 
         controls = ~ gdp + growth + trade + inflation + iso3c,
         data = bh_4y,
         number = 11)

# BEA CAPITAL EXPENDITURES
estimate(dv = 'bea_capitalexpenditures', 
         controls = ~ gdp + growth + trade + inflation + iso3c,
         number = 13,
         data = bh_4y)

# BEA TOTAL ASSETS 
estimate(dv = 'bea_totalassets', 
         controls = ~ gdp + growth + trade + inflation + iso3c,
         number = 15,
         data = bh_4y)

# REGION DUMMIES
estimate(dv = 'unctad_fdi_flow', 
         controls = ~ gdp + growth + trade + inflation + iso3c + region,
         number = 17,
         data = bh_4y)

# PPE COUNTRY-YEAR
estimate(dv = 'bea_ppe', 
         controls = ~ gdp + growth + trade + inflation + iso3c,
         number = 19,
         data = bh_1y)

# MULTIPLE IMPUTATION WITH LITERACY + TELEPHONE + URBAN
set.seed(1024)
dat_mi = amelia(bh_4y, m=5, cs='iso3c', ts='year', idvars='region', p2s=0)
apolitical = unctad_fdi_flow ~ gdp + growth + trade + inflation + gdppc + urban + literacy + telephone + iso3c
political = unctad_fdi_flow ~ gdp + growth + trade + inflation + gdppc + urban + literacy + telephone + political + iso3c
results = list()
results[['baseline']] = zelig(formula=apolitical, model='ls', data=dat_mi)
vars = grep('prs_', colnames(bh_4y), value=TRUE)
for(v in vars){
    for(i in 1:length(dat_mi$imputations)){ 
        dat_mi$imputations[[i]]$political = dat_mi$imputations[[i]][, v]
    }
    results[[v]] = zelig(formula=political, model='ls', data=dat_mi)
}

r2 = function(model){
    x = model$zelig.out$z.out
    x = sapply(x, function(i) summary(i)$r.squared)
    x = round(mean(x), 3)
    return(x)
}
dict = c('bureau'='Bureaucracy', 'corrupt'='Corruption',
         'demoacct'='Accountability', 'ethnic'='Ethnic Conf.',
         'extconf'='External Conf.', 'intconf'='Internal Conf.',
         'invprof'='Investment Profile', 'law'='Law and Order',
         'military'='Military', 'religious'='Religious')
names(dict) = paste('prs_', names(dict), sep='')
dict['baseline'] = 'Baseline'
r2 = sapply(results, r2)
names(r2) = names(results)
names(r2) = dict[names(r2)]
r2 = r2[!is.na(names(r2))]
r2 = data.frame('Variable'=names(r2), 'R2'=r2)
write.csv(r2, file='tables/table_8.csv')
